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Claims 

1. A method for seismic exploration, comprising actions of: 

(a) obtaining a set of seismic data traces representing a three-dimensional 
volume of seismic data samples; 

(b) dividing the three-dimensional volume into a plurality of smaller 
subvolumes; 

(c) selecting a discrete set of dip values and azimuth values; 

(d) dividing the three-dimensional volume into a plurality of parallelipipeds, 
each of the parallelpipeds being tilted by one of the selected dip values and 
rotated by one of the selected azimuth values; 

(e) halving each parallilepiped to obtain two half-parallelpipeds, wherein the 
two half-parallelpipeds have contain an equal number of samples and there 
exists a one-to-one relationship between corresponding samples in the two 
half-parallelpipeds; 

(f) enumerating the samples in each of the two half-parallelpipeds so as to 
obtain two vectors, such that corresponding samples in the two half- 
parallelpipeds have corresponding indices in the two vectors; 

(g) calculating a three-dimensional edge detection measure from the two 
vectors; 

(h) associating the computed edge detection measure to the parallilepiped 
center point to obtain a first subresult; 

(i) applying a constrast enhancement measure to the first subresult to obtain a 
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second subresult; 

(j) filtering the second subresult to obtain a set of third subresults by convolving 

the second subresult with a directional filter kernel that is tilted and rotated in 

accordance with a set of dip and azimuth values that correspond to the 

selected dip values and azimuth values of the computational analysis 

parallilepiped associated with the second subresult; 

(k) selecting a maximum filtered value from the set of third subresults; 

(I) applying a three-dimensional skeletonization algorithm to the maximum 

filtered value to generate a skeleton representing a fault surface; 

(m) executing action (I) with respect to the maximum filtered value for each of 

the plurality of parallelpipeds to obtain a plurality of distinct skeletons; and 

(n) labeling each of the plurality of distinct skeletons as a separate geologic 

feature. 

2. The method of claim 1_ , wherein the method comprises an additional action of: 
(b1) computing, for at least one of the subvolumes, a local seismic data 
orientation and an associated statistical uncertainty, and further, 

wherein each of the plurality of parallelpipeds is associated with a subvolume 
from the plurality of subvolumes and each of the plurality of parallelpipeds has 
a top surface and a bottom surface that are parallel to the local seismic data 
orientation for the subvolume associated with that parallelpiped. 

3. The method of claim 2, where action (b1 ) further includes: 
(b1a) applying a smoothing function to a subvolume; 
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(bib) computing, in three dimensions, a local orientation of seismic data within 
the subvolume; 

(b1c) averaging, over a surrounding three-dimensional space, the local 
orientation of seismic data within the subvolume to obtain a stable orientation 
estimate; and 

(bid) computing a measure of statistical uncertainty of the orientation estimate. 

4. The method of claim 2, where action (b1) further includes: 

(b1a) transforming the seismic data of a subvolume so as to map Cartesian 
coordinates of the subvolume into complex values in a complex plane; 
(bib) computing a three-dimensional gradient of the subvolume; 
(b1c) convolving the three-dimensional gradient with a three-dimensional 
Gaussian filter to obtain a convolution result; 

(bid) extracting a local orientation estimate of the subvolume as the argument 

(in the sense of a complex number) of the convolution result; and 

(b1e) computing an uncertainty measure (p) of the local orientation estimate as 

.where Ao and /h are eigenvalues of a scatter matrix for the subvolume's local 
spectrum. 

5. The method of claim 1^ wherein the selected dip values are selected from an 
interval extending from about -45 degrees to about 45 degrees and with a dip 
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increment selected from an interval extending from about 5 degree to about 10 
degrees. 

6. The method of claim 1_ f wherein the selected azimuth values are selected from 
an interval extending from about 0 degrees to about 180 degrees and wherein 
the selected azimuth values are selected using an azimuth increment selected 
from an interval extending from about 5 degrees to about 30 degrees. 

7. The method of claim 1_, wherein at least one of the parallelpipeds is tilted by a 
selected dip value and rotated by a selected azimuth angle to form a 
parallelpiped having integer dimensions of the form U x (2Z.2+1) * N, with N 
representing a number of samples measured with respect to a vertical 
dimension of the parallelpiped, with the selected dip value measured with 
respect to the vertical direction, and the azimuth measured with respect to a 
direction of north in the seismic data. 

8. The method of claim 7, wherein each parallelpiped has a minimum dimension 
Z.1 of 6, a minimum dimension 2Z.2+1 of 5 samples and a minimum vertical 
dimension N of 41 samples. 

9. The method of claim 7, wherein the set of selected azimuth values is 
exclusively dependent on U. 

10. The method of claim 9, where the set of selected azimuth values has a 
cardinality of at least 8. 

1 1 . The method of claim 1^ wherein the edge detection measure for a given 
parallelpiped having a central analysis point A is at least a function of 
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||v u (^-v 2t ,(r,Ht 
and 

||vuMIHI v mMI 

where v-u and V2,a are the two vectors generated from halving the 
parallelpiped, yand (pare the parallelpiped's dip, and azimuth, respectively, 
and p denotes a choice of vector space norm for a vector space in which vi 
and V2 are defined. 

12. The method of claim JM , wherein the three dimensional edge detection 
measure is a function of 

|vm(^)-v 2 ,.M1, 

13. The method of claim 12, wherein values of 

||v u (r^)-v 2> ,(r^)|| p 

and 

|vmM|,H|v m MI 

are computed only once per each pair of vectors v-w and V2,a and dip/azimuth 
combination, are subsequently stored in memory so as to memoize the 
computed values such that the computed values need not be recomputed for 
that pair of vectors and dip/azimuth combination. 

14. The method of claim 13, wherein the memoized computed values from a 
window of preceding edge detection measure computations are utilized, 
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through dynamic programming, to compute subsequent edge detection 
measures. 

15. The method of claim _1 , wherein the contrast enhancement measure is a 
function of a filter that is applied by convolving the first subresult with a rotated 
form of a "Mexican hat function" of the form 

f(n) = C(l-n 2 )e-" 2/ \ 

(where e is the base of the natural logarithm function) for values of n taken 
from an interval extending from about -4.5 to 4.5, the filter containing an odd 
number of filter coefficients and having its main direction perpendicular to slabs 
of the first subresult's corresponding parallelpiped, wherein convolving the first 
subresult with the rotated form of the "Mexican hat function" results in a 
"Mexican hat function convolution result." 

16. The method of claim 15, wherein the contrast-enhanced three-dimensional 
edge detection measure is computed as the maximum of the "Mexican hat 
function convolution result" and zero. 

17. The method of claim 1_, where portions of fault surfaces that are approximately 
aligned with a parallelpiped are extracted. 

18. The method of claim 17, wherein the directional filter kernel represents a three- 
dimensional finite impulse response filter having a set of coefficients 
characterized by a primary axis, a secondary axis, and a tertiary axis, wherein 
the secondary and tertiary axes are significantly shorter in length than the 
primary axis. 
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19. The method of claim 18, wherein the three-dimensional filter is rotated by an 
angle y+ a with respect to a time axis and rotating it by an angle 0 around an 
in-line axis, wherein ^represents a being a tolerance of dip increment. 

20. The method of claim 19, wherein the three-dimensional filter has a number of 
coefficients in the primary direction that is at least equal to a number of 
samples in a corresponding direction of the parallilepiped. 

21. The method of claim 20, wherein the three-dimensional filter is a long three- 
dimensional Hanning window with very small number of coefficients along the 
secondary and tertiary axes in relation to the number of coefficients along the 
primary axis. 

22. The method of claim 21^, wherein obtaining the set of third subresults includes 
applying, to a result of the three-dimensional filter, a threshold function defined 
by a user-determined threshold. 

23. The method of claim 22, further comprising filtering the set of third subresults 
to produce a three-dimensional volume output in which portions of fault 
surfaces parallel to the selected dip and azimuth values are distinguished. 

24. The method of claim 23, further comprising selecting the maximum 
directionally filtered output value over the filtered set of third subresults. 

25. The method of claim 24, wherein the three-dimensional skeletonization 
algorithm is applied to the maximum directionally filtered output value. 

26. The method of claim 24, wherein action (m) further includes actions of: 
(ml) computing a local orientation volume for the set of third results, the 
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computation of the local orientation volume being accomplished by computing 
local gradient vectors at each volume point using a 2-point stencil in each 
direction and finding a principal component of the local gradient vectors using 
singular value decomposition; 

(m2) finding the location xoi a largest-valued point in an output of action (I) 
that has not already been marked as masked point and ending execution of 
action(m) if largest-valued point's value is smaller than a threshold fault value; 
(m3) selecting points that are less than a predetermined distance of M from x 
in a horizontal plane containing x, with N\ and that have a local orientation 
vector that is within a preselected orientation difference threshold; 
(m4) selecting points that lie in a plane perpendicular to the local orientation 
vector at the point x, that are less than a predetermined distance of M from x, 
and that have a local orientation vector that is within a preselected orientation 
difference threshold; 

(m5) repeating actions (m3) and (m4), considering one iteration for all of the 
points selected in actions (m3) and (m4) until the total number of iterations 
exceeds a preselected number of iterations; 

(m6) marking output volume points around points of maximum fault volume 
value as masked points in response to a determination that a total number of 
points selected is less than a preselected number of minimum points; 
(m7) skeletonizing fault surfaces composed of the selected points from actions 
(m3) and (m4) in response to a determination that the total number of selected 
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points is greater than a preselected number of points; 
(m8) zeroing selected fault output volume points in response to a 
determination that a number of selected points is smaller than a preselected 
number of points; 

(m9) in response to a determination that a number of points selected during a 
immediately preceding iteration of skeletonization is greater than a preselected 
number of points, labelling the skeletonized fault surfaces with labels that are 
unique to each fault surface. 

27. The method of claim 26, wherein skeletonization includes: 

selecting a point "/ in accordance with either action (m3) or action (m4); 
searching along a line that passes through j/and is parallel to the local 
orientation vector and finding a location of a maximum fault output volume 
along that line that is within Nz grid points of y, 

marking in a new output volume the location of the maximum fault output 
volume value with a non-zero value; and 

removing j/from the set of points selected in actions (m3) and (m4) 

28. The method of claim 25, wherein the skeletonized filtered volume is stored in 
magnetic media for subsequent use. 

29. The method of claim 29, wherein descriptions of individual fault surfaces are 
stored in magnetic media for subsequent use. 

30. An apparatus for processing and analyzing seismic trace data, comprising 
means for: 
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(a) obtaining a set of seismic data traces representing a three-dimensional 
volume of seismic data samples; 

(b) dividing the three-dimensional volume into a plurality of smaller 
subvolumes; 

(c) selecting a discrete set of dip values and azimuth values; 

(d) dividing the three-dimensional volume into a plurality of parallelipipeds, 
each of the parallelpipeds being tilted by one of the selected dip values and 
rotated by one of the selected azimuth values; 

(e) halving each parallilepiped to obtain two half-parallelpipeds, wherein the 
two half-parallelpipeds have contain an equal number of samples and there 
exists a one-to-one relationship between corresponding samples in the two 
half-parallelpipeds; 

(f) enumerating the samples in each of the two half-parallelpipeds so as to 
obtain two vectors, such that corresponding samples in the two half- 
parallelpipeds have corresponding indices in the two vectors; 

(g) calculating a three-dimensional edge detection measure from the two 
vectors; 

(h) associating the computed edge detection measure to the parallilepiped 
center point to obtain a first subresult; 

(i) applying a constrast enhancement measure to the first subresult to obtain a 
second subresult; 

0) filtering the second subresult to obtain a set of third subresults by convolving 
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the second subresult with a directional filter kernel that is tilted and rotated in 

accordance with a set of dip and azimuth values that correspond to the 

selected dip values and azimuth values of the computational analysis 

parallilepiped associated with the second subresult; 

(k) selecting a maximum filtered value from the set of third subresults; 

(I) applying a three-dimensional skeletonization algorithm to the maximum 

filtered value to generate a skeleton representing a fault surface; 

(m) executing action (I) with respect to the maximum filtered value for each of 

the plurality of parallelpipeds to obtain a plurality of distinct skeletons; and 

(n) labeling each of the plurality of distinct skeletons as a separate geologic 

feature. 



